// Copyright 2007 by Rory Plaire (codekaizen@gmail.com)

using System;
using System.Collections.Generic;
using System.Text;

namespace Novas.Net
{
	public class CatalogEntry
	{
		private string _catalogId;
		private string _starName;
		private long _starNumber;
		private double _rightAscension;
		private double _declination;
		private double _properMotionRightAscension;
		private double _properMotionDeclination;
		private double _parallax;
		private double _radialVelocity;

		private CatalogEntry() { }

		public CatalogEntry(string catalogId, string starName,
			long starNumber, double rightAscension, double declination,
			double properMotionRightAscension, double properMotionDeclination,
			double parallax, double radialVelocity)
		{
			if (String.IsNullOrEmpty(catalogId)) throw new ArgumentNullException("catalogId");
			if (String.IsNullOrEmpty(starName)) throw new ArgumentNullException("starName");

			CatalogId = catalogId;
			StarName = starName;
			StarNumber = starNumber;
			RightAscension = rightAscension;
			Declination = declination;
			ProperMotionRightAscension = properMotionRightAscension;
			ProperMotionDeclination = properMotionDeclination;
			Parallax = parallax;
			RadialVelocity = radialVelocity;
		}

		public string CatalogId
		{
			get { return _catalogId; }
			private set { _catalogId = value; }
		}

		public string StarName
		{
			get { return _starName; }
			private set { _starName = value; }
		}

		public long StarNumber
		{
			get { return _starNumber; }
			private set { _starNumber = value; }
		}

		public double RightAscension
		{
			get { return _rightAscension; }
			private set { _rightAscension = value; }
		}

		public double Declination
		{
			get { return _declination; }
			private set { _declination = value; }
		}

		public double ProperMotionRightAscension
		{
			get { return _properMotionRightAscension; }
			private set { _properMotionRightAscension = value; }
		}

		public double ProperMotionDeclination
		{
			get { return _properMotionDeclination; }
			private set { _properMotionDeclination = value; }
		}

		public double Parallax
		{
			get { return _parallax; }
			private set { _parallax = value; }
		}

		public double RadialVelocity
		{
			get { return _radialVelocity; }
			private set { _radialVelocity = value; }
		}

		public CatalogEntry TransformEpoch(double currentDate, double newDate, string newCatalogId)
		{
			return Transform(CatalogEntryTransformType.ChangeEpoch, this, currentDate, newDate, newCatalogId);
		}

		public CatalogEntry TransformEquatorAndEquinox(double currentDate, double newDate, string newCatalogId)
		{
			return Transform(CatalogEntryTransformType.ChangeEquatorAndEquinox, this, currentDate, newDate, newCatalogId);
		}

		public CatalogEntry TransformEpochEquatorAndEquinox(double currentDate, double newDate, string newCatalogId)
		{
			return Transform(CatalogEntryTransformType.All, this, currentDate, newDate, newCatalogId);
		}

		public static CatalogEntry Transform(CatalogEntryTransformType transformType, CatalogEntry inputEntry, double inputDate, double newDate, string newCatalogId)
		{
			short j;

			double jd_incat, jd_newcat, paralx, dist, r, d, cra, sra, cdc, sdc,
			   term1, pmr, pmd, rvl, xyproj;

			double[] pos1 = new double[3];
			double[] vel1 = new double[3];
			double[] pos2 = new double[3];
			double[] vel2 = new double[3];

			/*
			   If necessary, compute Julian dates.
			*/

			if (inputDate < 10000.0)
			{
				jd_incat = Constants.Epoch + (inputDate - 2000.0) * 365.25;
			}
			else
			{
				jd_incat = inputDate;
			}

			if (newDate < 10000.0)
			{
				jd_newcat = Constants.Epoch + (newDate - 2000.0) * 365.25;
			}
			else
			{
				jd_newcat = newDate;
			}

			/*
			   Convert input angular components to vectors

			   If parallax is unknown, undetermined, or zero, set it to 1.0e-7 
			   second of arc, corresponding to a distance of 10 megaparsecs.
			*/

			paralx = inputEntry.Parallax;

			if (paralx <= 0.0)
			{
				paralx = 1.0e-7;
			}

			/*
			   Convert right ascension, declination, and parallax to position 
			   vector in equatorial system with units of AU.
			*/

			dist = Constants.RadiansToSeconds / paralx;
			r = inputEntry.RightAscension * 54000.0 / Constants.RadiansToSeconds;
			d = inputEntry.Declination * 3600.0 / Constants.RadiansToSeconds;
			cra = Math.Cos(r);
			sra = Math.Sin(r);
			cdc = Math.Cos(d);
			sdc = Math.Sin(d);
			pos1[0] = dist * cdc * cra;
			pos1[1] = dist * cdc * sra;
			pos1[2] = dist * sdc;

			/*
			   Convert proper motion and radial velocity to orthogonal components 
			   of motion, in spherical polar system at star's original position, 
			   with units of AU/day.
			*/

			term1 = paralx * 36525.0;
			pmr = inputEntry.ProperMotionRightAscension * 15.0 * cdc / term1;
			pmd = inputEntry.ProperMotionDeclination / term1;
			rvl = inputEntry.RadialVelocity * 86400.0 / Constants.KmPerAU;

			/*
			   Transform motion vector to equatorial system.
			*/

			vel1[0] = -pmr * sra - pmd * sdc * cra + rvl * cdc * cra;
			vel1[1] = pmr * cra - pmd * sdc * sra + rvl * cdc * sra;
			vel1[2] = pmd * cdc + rvl * sdc;

			/*
			   Update star's position vector for space motion (only if 'transformType' == CatalogEntryTransformType.Epoch
			   or 'transformType' == CatalogEntryTransformType.All).
			*/

			if ((int)(transformType & CatalogEntryTransformType.ChangeEpoch) != 0)
			{
				for (j = 0; j < 3; j++)
				{
					pos2[j] = pos1[j] + vel1[j] * (jd_newcat - jd_incat);
					vel2[j] = vel1[j];
				}
			}
			else
			{
				for (j = 0; j < 3; j++)
				{
					pos2[j] = pos1[j];
					vel2[j] = vel1[j];
				}
			}

			/*
			   Precess position and velocity vectors (only if 'transformType' == CatalogEntryTransformType.ChangeEquatorAndEquinox
			   or 'transformType' == CatalogEntryTransformType.All).
			*/

			if ((int)(transformType & CatalogEntryTransformType.ChangeEquatorAndEquinox) != 0)
			{
				for (j = 0; j < 3; j++)
				{
					pos1[j] = pos2[j];
					vel1[j] = vel2[j];
				}

				pos2 = Computations.Precession(jd_incat, pos1, jd_newcat);
				vel2 = Computations.Precession(jd_incat, vel1, jd_newcat);
			}

			/*
			   Convert vectors back to angular components for output.

			   From updated position vector, obtain star's new position expressed 
			   as angular quantities.
			*/

			xyproj = Math.Sqrt(pos2[0] * pos2[0] + pos2[1] * pos2[1]);
			r = Math.Atan2(pos2[1], pos2[0]);
			d = Math.Atan2(pos2[2], xyproj);

			CatalogEntry newEntry = new CatalogEntry();

			newEntry.RightAscension = r * Constants.RadiansToSeconds / 54000.0;
			newEntry.Declination = d * Constants.RadiansToSeconds / 3600.0;

			if (newEntry.RightAscension < 0.0)
			{
				newEntry.RightAscension += 24.0;
			}

			dist = Math.Sqrt(pos2[0] * pos2[0] + pos2[1] * pos2[1] + pos2[2] * pos2[2]);

			paralx = Constants.RadiansToSeconds / dist;
			newEntry.Parallax = paralx;

			/*
			   Transform motion vector back to spherical polar system at star's 
			   new position.
			*/

			cra = Math.Cos(r);
			sra = Math.Sin(r);
			cdc = Math.Cos(d);
			sdc = Math.Sin(d);
			pmr = -vel2[0] * sra + vel2[1] * cra;
			pmd = -vel2[0] * cra * sdc - vel2[1] * sra * sdc + vel2[2] * cdc;
			rvl = vel2[0] * cra * cdc + vel2[1] * sra * cdc + vel2[2] * sdc;

			/*
			   Convert components of motion to from AU/day to normal catalog units.
			*/

			newEntry.ProperMotionRightAscension = pmr * paralx * 36525.0 / (15.0 * cdc);
			newEntry.ProperMotionDeclination = pmd * paralx * 36525.0;
			newEntry.RadialVelocity = rvl * Constants.KmPerAU / 86400.0;

			/*
			  Take care of zero-parallax case.
			*/

			if (newEntry.Parallax <= 1.01e-7)
			{
				newEntry.Parallax = 0.0;
				newEntry.RadialVelocity = inputEntry.RadialVelocity;
			}

			/*
			   Set the catalog identification code for the transformed catalog
			   entry.
			*/

			newEntry.CatalogId = newCatalogId;

			/*
			   Copy unchanged quantities from the input catalog entry to the
			   transformed catalog entry.
			*/

			newEntry.StarName = inputEntry.StarName;
			newEntry.StarNumber = inputEntry.StarNumber;

			return newEntry;
		}
	}
}
